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A NEW LANCZOS-TYPE ALGORITHM FOR SYSTEMS 
OF LINEAR EQUATIONS 

MUHAMMAD FAROOQi AND ABDELLAH SALHR 


Abstract. Lanczos-type algorithms are efficient and easy to imple¬ 
ment. Unfortunately they breakdown frequently and well before con¬ 
vergence has been achieved. These algorithms are typically based on 
recurrence relations which involve formal orthogonal polynomials of 
low degree. In this paper, we consider a recurrence relation that has 
not been studied before and which involves a relatively higher degree 
polynomial. Interestingly, it leads to an algorithm that shows su¬ 
perior stability when compared to existing Lanczos-type algorithms. 
This new algorithm is derived and described. It is then compared 
to the best known algorithms of this type, namely A 5 /H 10 , As/^io, 
as well as Arnoldi’s algorithm, on a set of standard test problems. 
Numerical results are included. 

Key words : Lanczos algorithm; Arnoldi algorithm; Systems of Lin¬ 
ear Equations; Eormal Orthogonal Polynomials 

AMS SUBJECT : Primary 65F10. 


1. Introduction and background 

The Lanczos algorithm, [281 l29l E], is an iterative process that has 
been primarily designed to calcnlate the eigenvalnes of a matrix. How¬ 
ever, it has fonnd a wide application in the area of Systems of Lin¬ 
ear Equations (SLE’s) where it now is a well established solver. Its 
attraction resides in its efficiency as it only involves vector-to-vector 
and matrix-to-vector products. Moreover, in exact arithmetic, it con¬ 
verges to the exact solution in at most n steps, where n is the dimension 
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of the problem, [29]. While efficiency is its strong point, stability is 
not. Indeed, it is well known to breakdown as orthogonality of the so 
called Lanczos vectors, generated dnring the solution process, is lost. 
Efforts to avoid this breakdown led to a flurry of papers particularly 
from Brezinski and his team, [S] E] El El EDI [121 Ell E3], and others, 
|11EIE51EII[2Q1EIIE31|211[3IIE21I331[351EEIE9]. 

Several Lanczos-type algorithms have been designed and among them, 
the famous conjugate gradient algorithm of Hestenes and Stiefel, |25] . 
when the matrix is Hermitian and the bi-conjugate gradient algorithm of 
Fletcher, [ 22 ], and the algorithm of Arnoldi, mm, in the general case. 

Lanczos-type algorithms are commonly derived from recurrence rela¬ 
tions typically using Formal Orthogonal Polynomials (FOP’s) of low de¬ 
gree, [2H1 El ESI EZ]. Recurrence relations using relatively higher degree 
FOP’s have not been investigated. Here, we set out to design an algo¬ 
rithm that is based on such recurrence relations and FOP’s, and study 
its properties and, in particular, its stability. 

1.1. The Lanczos Process. Consider the SLE 

Hx = b, 

where A G b G R” and x G R". 

Let xq and y be two arbitrary vectors in such that y 7 ^ 0 
Lanczos method [29] consists in constructing a sequence of vectors 
defined as follows 

Xfc - Xo G Kk{A, ro) = span{ro, Hro,..., A’^~^ro), (2) 

Tfc = b - Hxfc-LLfc(H'^, y) = span{y, H'^y,..., (3) 

where A^ denotes the transpose of A. 

Equation (2) gives, 

Xfc - Xq = -oiro - Q! 2 Aro-afcH*’"^ro. (4) 

Multiplying both sides by A and adding and subtracting b on the left 
hand side gives 

rk = ro + oiro a 2 Aro H-h OfcH^'^ro. (5) 

If we set 

Pkix) = 1 -F aix + ... -t- akX, 
then we can write from (5) 


( 1 ) 

then 
Xfc G 


Ffe ^ Pt(A)ro. 


(6) 
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From (3), the orthogonality condition gives 

(yl'^V, rfc) = (y, = (y, A"Vk{A)vo) = 0, for i = 0,A; - 1. 

Thus, the coefficients ai,...,afc form a solution of SLE’s, 

ai(y, ^*+Vo) + ...+afc(y, =-(y, ^*ro), for z = 0,..., fc-1. (7) 

If the determinant of the above system is not zero then its solution exists 
and allows to obtain and r^. Obviously, in practice, solving the above 
system directly for increasing values of k is not feasible; k is the order 
of the iterate in the solution process. We shall see now how to solve this 
system for increasing values of k recursively, that is, if polynomials Pk can 
be computed recursively. Such computation is feasible, since polynomials 
Pk form a family of FOP’s which will briefly be explained below. 


1.2. Formal Orthogonal Polynomials. Dehne a linear functional c 
on the space of reel polynomials by 

c(a;*) = Ci for z = 0,1 ,... 


where 

Ci = (yl'^V, rfc) = (y, TTfc) for z = 0, 1,... 
Write the orthogonality condition as. 


c{x^Pk) = 0 for z = 0 ,..., fc — 1. 


( 8 ) 


The above condition shows that Pk is the polynomial of degree at most 
k and is a FOP with respect to the functional c, [S]. The normalization 
condition for these polynomials is Pk{0) = 1; exists and is unique if 
the following Hankel determinant 



Cl 

C2 • • 

Ck 

C2 

C3 • • 

Ck+l 

C/c 

C/c + 1 

• C2fc_i 


is not zero. In that case we can write Pk{x) as follows. 
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1 

X ■ ■ 

h 

X 

Co 

Cl • • 


^k—1 

Cfc 

* C2k-1 


Pk{x) = - 

Cl * * ‘ Cfc 


(9) 


Cfe • • • C2k-1 

where the denominator of this polynomial is the determinant of 
the system (7). We assume that V k, ^ 0 and therefore all the 
polynomials exist for all k. If for some k, = 0, then Pk does not 
exist and breakdown occurs in the algorithm, [a [ini dB da i. 


A Lanczos-type method consists in computing P^ recursively, then 
and dually x^. such that r^. = b — Ax^, without inverting matrix A. In 
exact arithmetic, this gives the solution of the system ( 1 ) in at most n 
steps, where n is the dimension of the SLE, mm- 


1.3. Notation and organization. The notation introduced by Baheux, 
in mm, for recurrence relations with three terms is adopted here. It puts 
recurrence relations involving FOP’s Pk{x) (the polynomials of degree at 
most k with regard to the linear functional c) and/or FOP’s pI^\x) (the 
polynomials of degree at most k with regard to linear functional 
M) into two groups: A, and Bj. Although relations A,, when they exist, 
rarely lead to Lanczos-type algorithms on their own (the exceptions being 
A 4 , [BE]; and Ai2, HZI, so far), relations Bj never lead to such algorithms 
for obvious reasons. It is the combination of recurrence relations Aj 
and Bj, denoted as Ai/Bj, when both exist, that lead to Lanczos-type 
algorithms. In the following we will refer to algorithms by the relation(s) 
that lead to them. Hence, there are, potentially, algorithms Aj and 
algorithms Ai/Bj, for some i = 1,2,... and some j = 1,2,... 

In this paper, a new algorithm based on a recurrence relation that has 
not been studied before, is derived. It is then compared to three other 
algorithms, one of which is the Arnoldi algorithm. 

The rest of the paper is organized as follows. In the next section, 
the Lanczos-type algorithm Ag/Hio, [2], and the estimation of the coeffi¬ 
cients of the recurrence relations Ag and Biq used to derive it are given. 
Section 3 is on the estimation of the coefficients of recurrence relation 
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Ai2, [E], used to derive the new algorithm of the same name. Section 
4 describes the test problems and reports numerical results. Section 5 is 
the conclusion and further work. 

2. Baheux algorithm A^/Biq 

The choice of algorithm A^/Biq^ for comparison with our own is dic¬ 
tated by the fact that this is the most robust of the algorithms considered 
in Eli on some of the problems considered here. So, outperforming this 
algorithm implies outperforming the rest of the algorithms considered 
therein. 

For completeness, we recall the relevant relations between adjacent 
FOPs that lead to A^/Biq and their coefficients estimates. These are As 
and Bio. The details of algorithm A^/Biq are given in [3]. 


2.1. Recurrence relation Ag. Relation Ag is 


Pk{x) = {AkX -f -h (CkX + Dk)Pk-i{x), 

(10) 

first investigated in El [3]. Its coefficients are estimated as 


Bk = 0, 

(11) 

Ck = 0, 

(12) 

Dk = 1, 

(13) 

and 


c{x'^-^ Pk-l{x)) 

Ak — 

c{x^pUi{x)) 

(14) 

As we know 


ic{x’‘Pk) = (A'^^y,Pfc(4l)ro) = (yfc,rfc), 

\c{x’^pI^^) = V)ro) = (yfc,Zfc), 

(15) 


with = A'^yf._^ and is defined in fl2l]) . Using flT^ . equation flTTj) 
becomes 


(yfc-i,rfc_i) ^ (yfc_i,rfc_i) 

(yfc,Zfc-i) (yfc-i,Tzfc_i)' 


(16) 
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2.2. Recurrence relation Biq. This relation, first investigated in [2l|3], 
is 

Pt\x) = (Alx + Bhd-I W + ClPtix), (17) 

Its coefficients are estimated as 


4 = 0, 

Clak = 1 , 

where Ofc is the coefficient of x*' in Pfc(x) dehned in flTOl) and 

oi _ 4c(yfc,rfc) 

— 7 r* 

c(yt.Zi-i) 


(18) 

(19) 


( 20 ) 


Equation 111 711 gives 

Zfc = Blzk-1 + Clrk. (21) 

2.3. Algorithm As/Biq. The pseudo-code of Ag/Pio, due to Baheux, 
PE], is as follows. 


Algorithm 1 Algorithm A^/Biq 

Choose xq and y such that y 7 ^ 0, and e arbitrarily small and positive. 
Set ro = b — Axq, 
zo = ro, 

To = Ti 

for /c = 0,1, 2,..., do 

A, — (yfcXfc) 

(ykA^k)’’ 

rfc+i = Tfc -|- Ak+iAzk, 

^fc+l Ak-\-lZk- 

if ||rfc+i|| > e, then 

yfc+i = A^y^, 

/^i 1 

k+l Ak+i ’ 

r>l _ _ ^k+1 (Yk+lAk+l) 

^+1 ~ (yfcPzfe) ’ 

Zfc+l = Pfc+iZfc -|- C^j^^Tk+l- 

else 

Stop; solution found. 

end if 
end for 
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3. A NEW Lanczos-type algorithm 

In the following, a new recurrence relation which leads to a new variant 
of the Lanczos algorithm is considered. 

3.1. Recurrence relation A 12 . Consider the following recurrence rela¬ 
tion, [mEg 

Pk(x) = Akl(x^ + BkX + Ck)Pk- 2 (x) + (BkX^ + EkX^ -h FkX + Gk)Pk- 3 (x)], 

(22) 

for k > 3, where Ak, Bk, Ck, Dk, Ek, Fk and Gk are constants to be deter¬ 
mined using the normalization condition Pk{0) = 1 and the orthogonality 
conditions c{x^Pk) = 0 , Vi = 0 ,..., A; — 1 , x* being a monic polynomial 
of exact degree i. To find these coefficients, we proceed as follows. Since 
Vfc, Pfc(O) = 1, equation (122]) gives 

1 = Afc [Gk + Gk] ■ 

Multiplying both sides of fl22]) by x* and then applying the linear func¬ 
tional c, we get 

c{x^Pk) = Afc{c(a;*+2Pfc-2) + Pfcc(P+^Pfc_ 2 ) Gkc{x^Pk- 2 ) + Dkcix^^^Pks) 

+EkC{x^~^‘^Pk-3) -|- FkC{x^~^^Pk-3) -|- GfcC(PPfc_3)}. 

(23) 

Equation (125]) is always true for i = 0,..., fc — 7. 

For i = k — 6, we have 

0 = Pfcc(a;^“^Pfc_ 3 ). 

Since c{x^~‘^Pk- 3 ) 7 ^ 0, we have 

Pfc = 0. 


For i = k — 5, we get 

0 = Ekc{x^-^Pk-3). 
But c{x^~'^Pk- 3 ) 7 ^ 0; therefore 

Ek = 0 . 


For i = k — A, (1251) gives 


Fk 


C{X^ ^Pfc-2) 
c{x^-^Pk-3) ' 


(24) 
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For i = k — 3, i = k — 2 and i 
respectively 

Bkc{x^-^Pk.2) + Gkc{x^-^Pk-^ 


fc — 1 we get the following equations 
---c{x^-^Pk-2)-Fkc{x^-^Pk-z). (25) 


BkC{x^-^Pk-2) + Ckc{x^-^Pk-2) + Gkc{x^-^Pk-^) 
= -c{x^Pk-2) - FkC{x^~^Pk-3), 

BkC{x^Pk-2) + GkC{x^~^Pk-2) + GkC{x^~^Pk-3) 
= -c{x’^^^Pk- 2 ) - FkC{x’^Pk-3)- 


(26) 

(27) 


Let Oil, cii 2 ) fli 3 ) 0215 022 ; 023 ; Osi; 032 ; sud CI 33 be the coefficients of Bk, 
Gk and Gk in equations fl25|l . fl26|) and fl27)) respectively and let bi, 62 and 
63 be the corresponding right sides of these equations. If represents 
the determinant of the coefficients matrix of the above mentioned system 
of equations then, we have 

ail = c{x’^~'^Pk-2), 012 = 0, Oi3 = c{x^-^Pk-3), 

021 = c(a;^“^Pfc_ 2 ), 022 = c{x^~‘^Pk-2), O23 = c{x^~‘^Pk-3), 

031 = c{x’^Pk-2), O 32 = c(a;''"^Pfc-2), O 33 = c(x^"^Pfc_ 3 ), 


bi = -c{x^-^Pk-2) - FkC{x^~‘^Pk- 3 ) = -O21 - Pfc 023 , 

62 = -c{x^Pk-2) - FkC{x^~^Pk-3) = -O31 - Pfc033, 

63 = -c(a;^+^Pfc_ 2 ) - Fkc{x^Pk-3) = -s - Fkt, where s = c(a;^+^Pfc_ 2 ) 
and t = c{x^Pk- 3 )- 

Therefore, equations (l25D . (I26|) and (HTh can be written as 


oiiPfc + OGk + cii3Gk = bi, (28) 

02lPfc + 022Cfc + Cl23Gk = b2, (29) 

031 -Bfc + 032 Cfc + a33Gk = b3. (30) 

To solve for Bk, Gk and Gk, Cramer’s rule requires 

Afc = 011(022033 — 032 O 23 ) + 013(021032 — 031022 )- 
If Ak 7 ^ 0, then 


Bk 


&l(022033 — O32O23) + 013(62032 — 63O22) 


Gk — 


Ak 

61 — OiiPfc 


Gk = 


Ol3 

62 — 02 lPfc — 023Gfc 


(31) 

(32) 


022 


(33) 
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and 

1 = [Cfc + • (34) 

With all the necessary coefficients now determined, the expression of 
the polynomials Pk{x) becomes 

Pkix) = Ak{{x^ + BkX + Ck)Pk-2{.x) + {FkX + Gk)Pk-3ix)}- (35) 

Let ns now nse the relation fl55]) to compnte Pk{x), necessary for the com- 
pntation of the residnal = b — = Pk{A)rQ and the corresponding 

vector Xfc. 

Assnme that Pk has exact degree k and the 3-term recnrrence rela¬ 
tionship (15^ holds. To move to the Krylov space, replace x hy A and 
mnltiply both side of (15B]) by tq to get, 

Pk{A)TQ = Ak[{A‘^+BkA+Ckk)Pk-2{A)ro + (FkA+Gkk)Pk-3{A)ro]. (36) 
Using eqnation (6), gives 

Tfc = Ak{{A^ + BkA -|- C'fc/)rfc_2 -|- {FkA + Gfci)rfc_3}. (37) 
And nsing = b — Axfc, gives 

Xfc = Ak{Gk^k-2 + Gk^k-3 - {Ark-2 + 5fcrfc-2 + Fk^k-s)}, (38) 

with Fk as in eqnation fl24D . Using (IT^ . Fk can be written as 

p _ (yfc-2) rfc- 2 ) 

^ (yfc_3, rfc- 3 ) ■ 

Condition flT5l) can be nsed eqnally to rewrite the expressions of an, 
throngh 033 , bi to 63 as follows. 

= (yfc-2Afc-2), ai2 = 0, ai3 = (y;i,_3, Ffc-s), 

021 = (yfc_l,rfc-2), 022 = (yA:-25 rfc- 2 ), O 23 = (yfc-25’’fc-s), 

031 = (yA„rfc_2), 032 = (yfc_i,rfc_2), 033 = (yfc_i,rfc_ 3 ), 

bi = —021 — Tfca23, ^2 = —031 — Fkttss, bs = —s — Fkt, 
where s = (y^+i, rfc_ 2 ) and t = (y^^,, r^.g). 

These parameters allow the explicit compntation of Bk, Gk, Gk, and Ak 
as is given by eqnations 31, 32, 33 and 34 respectively. Eqnations (l37)l 
and (l38ll define the new Lanczos-type algorithm. 
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Now, since all previous formulae are only valid for /c > 3, it is necessary 
to find the expressions of the polynomials of degrees 1 and 2. From (9), 
we can write 


Pi{x) 


1 X 

Co Cl 
Cl 


Pl{x) = 1 — —X, 

Cl 

ri = To - ^Aro, and xi = xq + ^ro, where c* = (y, yfVo). 
Using (9) again, we can write 


1 

X 


Co 

Cl 

C2 

Cl 

C2 

C 3 


Cl C2 
C2 Cs 


P2{x) 


1 - 


CqCs — C1C2 

- T^X + 

C1C3 - C2 


C0C2 Cl 2 

-2^ ’ 

C1C3 - C2 


r 2 = I’D — ct^ro + /d^^ro, and X 2 = xq + ctro — /d^ro, where ct = 
/d = and S = C 1 C 3 — C 2 . 


3.2. Algorithm A 12 . Putting together the various steps given in the 
above section, the new algorithm can be described as follows. 
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Algorithm 2 Algorithm 

Choose xo and y such that y 7 ^ 0, and choose e arbitrarily small and 
positive. 

Set ro = b - Axq, yo = y, P = ^ro, Pi = Ap, cq = (y, fq), 

Cl = (y, p), C 2 = (y, Pi), C 3 = (y, Ap^), 6 = cics - c\, 

2 

_ CQCS—C1C2 Q _ C 0 C 2 —Cl 

5 1 H s ^ 

ri = ro - ^p, xi = Xo + ^ro, 

r2 = ro - op + /3pi, X 2 = xo + aro - /3p, 


Yi = y2 = Ys = ^^y2, k = 3. 
while I |rfc| I > e do 

Yfc+i = A^y;., = Arfc_i, ^3 = ^^k-2, 

®ii = (yfc-25rfc-2), ai3 = (yfc_3,rfc_3), 021 = (yfc_i,rfc_2), 022 = an, 
a 23 = (yfc -25 i’fc- 3 ), a 3 i = (yfc,rfc_2), 032 = 021, 033 = (yfc_;L, rfc_3), 
s = (yfc+i,rfc_2), t = (yfc,rfc_3), Fk = 

bi = —021 — o,23Fk, 62 = —031 — a^sFk, 63 = —s — tF^, 

= 011(022033 — 032023) + 013(021032 — 031022), 

D _ bi{a22a,33—a32a23)+ai3{b2a32—b3a22) 

-, 

Q _ bi—aiiilfc 

^ ai3 ’ 

_ fe2—021 life—023Gfc 


022 


Ck = 

A, — ^ 

Cfe+Gfe ’ 

Tfe = Afc{q2 + i^fcqi + CkTCk-2 + -^fcQs + Gfcrfc- 3 }, 

Xfc = Afc{CfcXfc_2 + GkXk-3 - (qi + BkYk -2 + FkTk-s)}, 


k = k + 1- 

end while 

Stop; solution found. 


4. Numerical results 

Ai 2 , [E], the algorithm described in the above section, has been tested 
against algorithms A^/Biq and A^/Biq, the best algorithms according to 
[ 3 li, as well as against the established Arnoldi algorithm, [DIM]. 


4 . 1 . Test problems I. The test problems considered here arise in the 5- 
point discretisation of the operator —on a rectangular region. 
Comparative results on instances of the following problem ranging from 
dimension 10 to 100 for parameter 6 taking values 0.0 and 0.2 respectively, 
are recorded in Table 1 and Table 2. 
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p -/ ■■■ 

\ -I B -I 


A = 


0 \ 


, with B 


/ 4 a 

1/9 4 a 


0 \ 


; -1 B -1 

\ 0 . -J B 

and ct = —1 + and (3 = —1 — 5. 
X = (1,1,, 1 )'^, is the solution 


J 


I : (3 A a \ 

V 0 ■■■ f3 A ) 


The right-hand side is b = Tx, where 
of the system. The dimension of B is 


10 . 


Table 1. Experimental results for problems when <5 = 0 



1 Arnoldi 

1 4 I 5 /Bio 

1 ^s/-Bio 

1 ^12 1 

n 

\\rk\\ 

t(sec) 

Ikfcll 

t(sec) 


t(sec) 

Ihd 

t(sec) 

10 

1.2514E-1U 

0.001450 

2.2940E-i‘' 

0.001892 

1.7704E-i‘' 

0.002770 

4.9623E-i‘* 

0.002819 

20 

1.7733E-11 

0.002207 

2.5256E-i'‘ 

0.001842 

1.7489E-13 

0.002654 

1.7536E-13 

0.002904 

Kfll 

1.2990E-i'‘ 

0.003602 

3.9026E-°'^ 

0.002220 

4.9472E-°'^ 

0.003179 

5.4705E-°® 

0.003370 

40 

3.5434E-11 

0.006071 

1.4770E-i° 

0.002416 

8.4658E-i° 

0.003095 

1.4776E-'>» 

0.003526 

50 

6.1827E'°® 

0.008870 

1.9959E-“® 

0.002962 

1.3598E-“® 

0.003696 

4.7994E'°® 

0.003980 

60 

2.9843E“i'‘ 

0.012282 

9.1910E-°® 

0.003001 

3.7470E-0® 

0.003776 

5.0010E-'>« 

0.004354 

70 

4.2642E-13 

0.017151 

4.9035E-“® 

0.003622 

4.2579E-“® 

0.004194 

1.3781E^°® 

0.005658 

80 

5.0951E-°® 

0.021938 

4.4311E-°® 

0.004498 

7.7199E-°^ 

0.005504 

7.5581E-°® 

0.005271 

90 

9.6960E-13 

0.029083 

NaN 


8.5560E-°® 

0.007900 

3.7301E-'>s 

0.006541 

100 

1.1397E^i® 

0.036462 

1.1889E-°® 

0.003849 

3.1695E-°® 

0.004499 

8.9530E-°'^ 

0.005084 


Table 2. Experimental results for problems when 5 = 0.2 



1 Arnoldi 

1 7 I 5 / Bio 

1 ^s/Bio 

1 ^12 1 

n 

\\rk\\ 

t(sec) 

Ikfcll 

t(sec) 

Ihd 

t(sec) 

Ihd 

t(sec) 

10 

2.3499E-1-’ 

0.001377 

5.2347E-^4 

0.002339 

5.2347^-"^ 

0.002948 

5.2347E-''’'‘ 

0.003231 

20 

5.6622E““ 

0.002149 

4.1778E-11 

0.001842 

5.8526E^“ 

0.003090 

6.3915E-i“ 

0.003372 

30 

6.8771E~i® 

0.003573 

8.9881E-°'‘ 

0.002220 

8.9880^-°'* 

0.003580 

8.9880E-“ 

0.003847 

40 

1.8106E-i“ 

0.006137 

8.7583E-°'‘ 

0.002830 

9.3988^-°'* 

0.003620 

9.1261E-“'‘ 

0.003977 

50 

3.5345E-°® 

0.008552 

6.2669E-°'‘ 

0.003360 

5.7269^-°'* 

0.004055 

2.5040E-“'‘ 

0.004964 

60 

2.8757E-13 

0.012544 

6.3670E-“'‘ 

0.003877 

8.4915E-°'‘ 

0.004885 

7.3489E-“'‘ 

0.005345 

70 

4.2552E“13 

0.017352 

8.5670E~°'‘ 

0.003902 

7.0703E-°'‘ 

0.006158 

9.9086E-“ 

0.005052 

80 

1.7785E-“'‘ 

0.021629 

NaN 


NaN 


6.5602E-“'‘ 

0.012131 

90 

1.4837E-“'‘ 

0.029332 

NaN 


7.5451E-°'‘ 

0.011230 

9.5294E-“'‘ 

0.011842 

100 

5.8942E^i3 

0.037067 

NaN 


NaN 


9.9710E-“'‘ 

0.018899 


4.2. Test problems II. The coefficient matrix here is taken as the 
Hilbert matrix, i.e. A = hilb{n), where hilb{n) is a Matlab function, 
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n being the dimension of A. The right-hand side b and the solution x, 
are dehned in the same way as in test problems I. The Hilbert matrix 
is notoriously ill-conditioned. Ill-conditioned systems of linear equations 
are notoriously difficult to solve to any useful accuracy, dSlEZlISl] 

Table 3. Experimental results when H is a Hilbert matrix. 



Arnoldi 

7I5/ Biq 

Biq 

^12 

n 

Ihdl 

t(sec) 

IKII 

t(sec) 

IKII 

t(sec) 

Ihdl 

t(sec) 

10 

UTTole^ 

0.002364 

4.4809E-'^“ 

0.001945 

4.4809E-““ 

0.002834 

4.4809^-*^“ 

0.002706 

20 

2.2288E"15 

0.003285 

9.6002E^'’5 

0.001772 

9.6002E^'’5 

0.002839 

9.6002E"°5 

0.002787 

30 

3.8953E-15 

0.004148 

1.3341E-°5 

0.001929 

1.3376E-°® 

0.003050 

1.3340E-°5 

0.002791 

40 

3.2251E~i'‘ 

0.005521 

3.8228E~°5 

0.001880 

3.8276E“® 

0.003327 

3.8229E~°5 

0.002772 

50 

2.8673E~i® 

0.007068 

7.9457E~°5 

0.001997 

7.9463E“°® 

0.003475 

7.9457E~°5 

0.003334 


All algorithms have been implemented in Matlab version 7.8.0 and 
run on a PC, under Microsoft Windows XP Professional Operating Sys¬ 
tem, with 3.2GB RAM, and 2.40 GHz Intel(R) Core(TM) 2 CPU 6600. 
The problems are solved as dense problems, i.e. no sparsity has been 
exploited. The results point to the Arnoldi algorithm being the most ro¬ 
bust overall, but also the slowest overall. Algorithms A^/Biq and Ag/i^io 
are the fastest overall, but the least robust overall; in fact they have 
failed to solve some problems in high dimension due to breakdown, of 
course, which is endemic in Lanczos-type algorithms. Algorithm A 12 , 
like Arnoldi, solved all problems but faster and not as accurately. It is 
also more robust than A 5 /R 10 and Ag/Rio overall, but slower than both 
of them overall. Its lower speed compared to that of A 5 /R 10 and Ag/Rio 
is expected since the recurrence relation A 12 involves more coefficients 
than both recurrence relations A^/Biq and Ag/Rio. Note that on the 
Hilbert-type problems. Table 3, the algorithms could not cope with di¬ 
mensions higher than 50. Arnoldi is again the most stable overall and the 
slowest as the dimension grows. The other three algorithms have similar 
performances. 


5. Conclusion and further work 

The way Lanczos-type algorithms are derived using recurrence rela¬ 
tions involving POP’s means that many such algorithms can be created, 
each based on a different set of relations. The choice of recurrence re¬ 
lations to use is dictated by the degree of POP’s to be involved; high 
degrees mean a large number of coefficients have to be calculated in the 
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concerned Lanczos process. This, consequently, dictates the computa¬ 
tional complexity of the resulting Lanczos-type algorithm. However, it 
is well known, [221 ED], that computational complexity does not always 
imply efficiency, or indeed robustness. Moreover, robustness and accu¬ 
racy are often more important. It is therefore worthwhile to look beyond 
complexity issues sometimes, like we did here. 

In this paper we have shown that, indeed, there are recurrence rela¬ 
tions worth exploring since they lead to more robust algorithms. As a 
result, a new Lanczos-type algorithm, A 12 has been designed. The nu¬ 
merical performance of this algorithm is compared to that of two existing 
Lanczos-type algorithms, which were found to be the best among a num¬ 
ber of Lanczos-type algorithms, Eli, on the same set of problems as 
considered here. It is also compared to the well established Arnoldi algo¬ 
rithm. It is interesting to hnd that algorithm A 12 is overall more robust 
than A 5 /H 10 and Ag/Hio and faster than Arnoldi’s. This makes it occupy, 
at least on the set of problems used here and elsewhere, a happy medium 
position. It is therefore the ideal candidate for time-limited applications 
which do not require high accuracy. 
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